Laterally driven interfaces in the three-dimensional Ising lattice gas 
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O 

^o We study the steady state of a phase-separated driven Ising lattice gas in three dimensions using 

computer simulations with Kawasaki dynamics. An external force field F(z) acts in the x direction 

parallel to the interface, creating a lateral order parameter current j x (z) which varies with distance 

z from the interface. Above the roughening temperature, our data for 'shear-like' linear variation of 

F(z) are in agreement with the picture wherein shear acts as effective confinement in this system, 

thus supressing the interfacial capillary-wave fluctuations. We find sharper magnetisation profiles 

and reduced interfacial width as compared to equilibrium. Pair correlations are more suppressed in 

the vorticity direction y than in the driving direction; the opposite holds for the structure factor. 

Lateral transport of capillary waves occurs for those forms of F(z) for which the current j w (z) is an 

odd function of z, for example the shear-like drive, and a 'step-like' driving field. For a V-shaped 

driving force no such motion occurs, but capillary waves are suppressed more strongly than for 

the shear-like drive. These findings are in agreement with our previous simulation studies in two 

dimensions. Near and below the (equilibrium) roughening temperature the effective-confinement 

I picture ceases to work, but the lateral motion of the interface persists. 
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I. INTRODUCTION 



I 

Dimensionality d of space is a relevant parameter in condensed matter systems that exhibit large spatial fluctuations, 

e.g., thermal fluctuations of the order parameter near a second-order phase transition. Such critical fluctuations are 

correlated over large distances, which gives rise to singular behaviour characterized by a set of critical exponents 
i ^i whose values depend on d. Moreover, the very existence of the phase transition depends on d. Below the lower critical 

dimension fluctuations are strong enough to destroy the ordered phase, and hence there is no phase transition. In 

contrast, above the upper critical dimension fluctuations are no longer important, and the critical exponents become 
^ the same as in mean field theory. 

An important dependence on dimensionality also occurs in systems where two distinct phases coexist. In such 

systems thermal fluctuations can be correlated over a long distance within the interfacial region. Then the interfacial 
_i- correlation length £y (parallel to the interface) increases with increasing thickness w of the interface, £m oc w^, 

where < £ < 1 lj, and the interface is termed rough. The value of the roughness exponent £ again depends on 
"iJ dimensionality d. A statistical-mechanical description of interfacial degrees of freedom, capillary wave theory, was 

proposed by Buff, Lovett and Stillinger [31 [3] , who suggested to model the interface as a sharp divide between the 
,__! two phases, but one that would freely fluctuate. Capillary wave theory (CWT) predicts that £ = ^(3 — d) for d < 3, 

i.e., fluctuations can destabilise the interface in d < 3. Then the thickness w of the free interface is infinite in the 

thermodynamic limit. At the marginal dimension d — 3, the value £ = corresponds to a logarithmic divergence. For 
S^ d > 3, one has ( — but w is finite for all £y ; the interface is then said to be smooth. 

These results originate from the choice of the simplest, Gaussian form of fluctuations, i.e., the probability for a local 

departure (height) h(r) of the interface from the reference plane ft, = is given by the Boltzmann weight ex e - A ' H /' c B T ; 

where ksT is the thermal energy, and 

AH = J d d ~ l v f E [v^r)] 2 + V«t[ft(r)]) , (1) 

where r are coordinates parallel to the interface. T is the interfacial stiffness, which for a continuum fluid is simply 
the tension of a free interface [TJ|4]-[6]. If the external potential V cx t is quadratic in h(r), for example due to gravity, 



the equipartition theorem for quadratic degrees of freedom may be applied. The equilibrium interface pair correlation 
function for a translationally-invariant system is then found to be 

n( \ lh(n\u\\ ksT f d ^ lq e ' qr a\ 

C(r) = (h(O)h(r)) = —j (27r)d _ le _ 2 + g2 , (2) 

where £7 2 = r _1 (9 2 V xt / dh 2 ) , and the limit of infinite interface dimension L d ~ x , L — > oo has been taken. An upper 
cutoff on wave numbers |q| < tt/cl is always assumed, a is usually identified with some appropriate microscopic length 
in the interface region, e.g., the lattice spacing or the bulk correlation length £& 7\. The behaviour of C(r) is obviously 
d-dependent - as is the interface width w, defined by w 2 = C(0). 

For certain microscopic models, such as the Ising model, which is equivalent to the lattice gas model of a fluid, exact 
results are available H] . A description starting from a microscopic Hamiltonian is particularly useful because it 
accounts for both interfacial and bulk fluctuations. In contrast bulk degrees of freedom of coexisting phases are absent 
in CWT, in which one considers only interfacial configurations. On the other hand, in the "classical" theories for 
the interface based on order parameter (or density) profiles, for example the van der Waals theory [9J [TU] , interfacial 
fluctuations are absent and the structure of the interface is reduced to an inhomogeneity of the order parameter, i.e., 
it is of order of the bulk correlation length £&. 

In d = 2, Ising interfaces are rough for all temperatures T below the critical temperature T c , and w 2 ex L, as revealed 
by exact results for the local magnetisation profile on a scale large compared to £b [IT], in agreement with CWT. 
However, in d = 3 there is evidence [5l [T2l [13] for the existence of a finite roughening temperature < Tr < T c , below 
which the interface is smooth. This evidence is supported by Monte Carlo (MC) simulations [H] and by rigorous 
analysis of discrete random surface models, such as the solid-on-solid (SOS) and discrete Gaussian (DC) models 
[131 [TS]. These models approximate the interface of the Ising model at low temperatures, for example the SOS model 
can be regarded as a certain anisotropic-coupling limit of the Ising model where the SOS interface configurations are 
selected from the Ising configurations by the requirement that there are no overhangs or bubbles |8|. 

In the situations that we shall address, the interface is stabilised by the presence of two walls at spacing L±. This 
problem is interesting because the fluctuating interface will experience collisions with the constraining walls. The 
standard capillary wave model does not apply for this case; it has to be extended to take into account entropic 
repulsion from the walls. Confined Ising interfaces can be treated rigorously in d — 2. Results for the magnetisation 
profile |16j indicate that in two dimensions, in spite of the entropic repulsion at its extremities, the interface meanders 
back and forth between the walls so that w ex Lj_. This is markedly different from what is expected in three 
dimensions on the basis of the analysis of low-energy excitations in discrete random surface models |15] . An energy- 
versus-entropy argument leads to the conclusion that interface configurations which result in interference with the 
boundary are needle- like. A rigorous analysis gives w 2 ex L± and L± oc ln£n- These results support conjectures from 
the phenomenological effective interface Hamiltonian [TJ [H [T7] , and for d = 3 Ising interfaces they agree with MC 
simulation studies [18] , 

Dependence of the structure of equilibrium interfaces on the spatial dimensionality d has repercussions also on 
relaxation dynamics of the fluctuating interface. Relaxation dynamics of capillary waves have been recently studied 
using molecular dynamics simulations of simple liquids in d = 3 [6l [19] . For the liquid interface constrained between 
two walls at separation Lj_ , a pronounced enhancement of the relaxation time of capillary waves was found; the most 
affected are relaxation times of waves with the wave number q ~ L~J_ |19j . 

These results for equilibrium suggest that a non-trivial dependence of the interfacial structure and dynamics on 
dimensionality may persist to non-equilibrium situations. In this paper, we study this problem for fluctuating interfaces 
that are driven into a steady state by the action of an external field parallel to the plane of the interface. The problems 
of roughness, spatial and temporal correlations of laterally driven interfaces were first addressed in the lattice gas 
driven by a spatially uniform force field via MC simulations [501 HI]- These studies, which were of two-dimensional 
systems, found that the interface becomes less rough when drive is applied. The order-parameter (magnetisation) 
profiles become much sharper upon increasing the drive and the interfacial width is reduced as compared to the 
equilibrium value. By a suitable coarse graining of microscopic particle configurations, the local position (height) 
of the interface was defined, and the behaviour of the spatial interface height correlation function was studied. The 
results are consistent with a reduction of the interfacial correlation length £y. Moreover, the structure factor S(q), 
defined as the the Fourier transform of the height-height correlation function displays deviations from the equilibrium 
capillary wave dependence 1/q 2 as q — > 0; the data suggest a weaker singularity S(q) oc 1/g 067 , which implies a 
reduction of the roughness exponent £. In theoretical attempts to treat driven interfaces, one derives a dynamic 
equation for the interfacial degrees of freedom, starting from the time-dependent Landau-Ginzburg-type equation for 
the order parameter [23J[S3]. This approach leads to a non-local and non-linear equation for the interface height. A 
linear stability analysis of this equation for a spatially uniform drive parallel to the interface shows that temporal 
decay of fluctuations along the driving field is faster than that orthogonal to the driving field. However, predictions 



for the roughness of the interface do not agree with the simulation results in d = 2. 

More recent interest in the theoretical challenges of laterally driven interfaces originates from an experiment on 
colloidal gas-liquid interfaces subjected to shear flow. In the experiment of Derks et al. [21] a real-space visualisation of 
interfacial fluctuations revealed a reduced intcrfacial roughness when shear was applied. The width w and correlation 
function of the height of the interface were analysed by fitting to the equilibrium CWT results. The fitting parameters 
in the analytic CWT expressions were the correlation length £m along the interface measured in the flow direction, 
and the surface tension. The authors concluded that w was decreased but £n was increased. Recently, the problem of 
non-equilibrium fluctuations of a liquid-liquid interface under shear has been addressed theoretically [2 5) within the 
framework of fluctuating hydrodynamics. This approach leads to a mode-coupling equation for the interface height 
which was solved using a perturbation theory. Results for the interfacial width are in agreement with the experiment 
of Ref. [21], but the results for the interfacial correlation length in the flow direction are not. The theoretical height- 
height correlation function and the structure factor imply a decrease of the correlation length in the direction of flow. 
Interestingly, in the direction perpendicular to the flow (vorticity direction), the correlation length seems to increase. 
Similar conclusions have been obtained from molecular dynamics simulations [19]. 

Previously we have studied interfaces in the two-dimensional Ising strip driven by an external field that is applied 
parallel to the walls (and to the interface), and may vary in the direction perpendicular to the mean position of 
the interface, by using MC simulations with spin-exchange Kawasaki dynamics [26H28] . These studies were partially 
motivated by the need to understand sheared fluid interfaces. Because our results were obtained in d = 2, and because 
of the simplified character of our model and its dynamics, we were not in the position to attempt a direct comparison 
with experimental data. However, our results were in partial qualitative agreement with Ref. [21]. We found that the 
shear-like drive acts as an effective confinement on the system; a steady state is reached in which the magnetisation 
profile is the same as that in equilibrium, but with a rescaled length implying a reduction of the interfacial width. 
Pair correlation functions along the interface decay more rapidly with distance under drive than in equilibrium, and 
for cases of weak drive can be rescaled to the equilibrium result. Moreover, we find that interfacial transport can 
occur in an unexpected way parallel to the interface. The lateral flux of the order parameter at a planar interface 
induces lateral motion of the thermal capillary waves, provided that the flux is an odd function of distance from the 
interface. 

In the present paper we study the same model system using the same approach, but in three spatial dimensions. 
We wish to investigate to what extent our findings from 2d persist to higher dimensions. Moreover, dependence on 
dimensionality of various quantities characterising the structure and dynamics of the interface, as well as of transport 
properties in a driven state is interesting. One would like to know whether the different character of interfacial 
fluctuations, i.e., "wandering" in d — 2 and "spikes" [T3] in d = 3 has any repercussions. There are also new aspects 
that deserve to be studied: the behaviour of the system near the roughening transition, and the anisotropy effects 
introduced by the introduction of the vorticity direction. Last but not least, results in d = 3 permit a more direct 
comparison to experiment. 

The rest of this paper is organised as follows. In Section [n] we introduce the model and give details of the 
simulations. In Section [ill A[ we present and discuss results for the interfacial structure of the driven 3d Ising system, 
via the magnetisation profile, interface width, and correlation functions in real and Fourier space. In Section |IIIB[ 
we investigate the dynamics of the driven interface, showing results for the current and evidence for capillary wave 



transport. Finally we draw conclusions in Section IV 



II. MODEL AND SIMULATION DETAILS 

We consider a three-dimensional (3rf) Ising model on a simple cubic lattice. On lattice sites i sit spins <7j, which may 
take values ±1. In lattice gas language, the spin variables become particle occupation numbers Tj = (a + l)/2 = 0, 1, 
corresponding to the absence and presence of a particle at a site, respectively. The Hamiltonian for the system is 

H = -jy^a i (T j> (3) 

(id) 

where (i,j) indicates that the sum is over nearest neighbour sites i and j. The coupling constant J > 0, so that the 
interactions are ferromagnetic (attractive in lattice gas language). 

The lattice has dimensions L x x L y x L Zl with periodic boundary conditions applied in the x and y direc- 
tions. (All lengths are expressed in units of the lattice constant b = 1). Spin layers are located at heights 
z = — L ' 2 ~ 1 , — L ' 2 ~ 3 ...,— 1,2,..., I " 2 ~ 1 , for a total of L z layers. The interface is induced by walls of fixed spins 
a = +1 at the top (z — (L z + l)/2) and a — — 1 at the bottom (z — —(L z + l)/2) planes of the lattice; these boundary 
conditions energetically favour parallel alignment of the interface with the x-y plane, with the '+' phase in the upper 



half of the volume, z > 0. We focus on slab-like lattice geometries, with L x ,L y ^> L z and L x = L y = L, so that the 
system is confined between the two walls, and the scaling length scale for the interfacial width is L z = L± [4| IT5] . 

Time evolution of the system proceeds under Kawasaki 29J spin-exchange dynamics, which conserve the magneti- 
sation, or equivalently the number of lattice-gas particles, locally. Elementary simulation moves consist of exchanging 
the values of two randomly chosen nearest-neighbour spins with probability p; in the lattice gas, this corresponds to 
a particle moving to a neighbouring empty lattice site. The key features of these dynamics are their conservation 
of order parameter, and their locality. These attributes are desirable from the point of view of simulating particle 
movement (albeit crudely) and the competition between diffusive motion and driven motion. In the general case, an 
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FIG. 1: (Color online) Illustration of the force field, microscopic particle motion, and coarse-grained interface motion in the 
model system, (a) Shear-like driving field F x (z) = 72. The 'intruders' of one phase into the other move in the same direction 
in both the upper and lower halves of the system. As explained later, the interface displays transport, (b) V-shaped driving 
field: the intruders now move in opposite directions. No interfacial motion occurs; interfacial fluctuations decay and new ones 
appear as time passes. 



external force field F(z) — (F x (z),F y (z),0) acts on the system, driving in the x-y plane. This field alters the Monte 
Carlo acceptance rates, to produce a modified Metropolis rate, 



p = min{l, exp [-/3(AH + AW)}} . 



(4) 



Here, (3 = 1/k^T is the inverse temperature (the Boltzmann constant is set to unity), and AH is the change in 
internal energy from the proposed exchange. AW is the work done by or against the external force field; for AW = 0, 
the above rate reduces to the standard Metropolis one, which samples thermal equilibrium states. We are interested 
in the case of non-zero AW, when the system will reach a non-equilibrium steady state. The system is immersed in a 
heat bath at constant temperature T, into which the work done is dissipated. The driving field is related to the work 
term by 



AW = -JS-F(z)(a l -a J )/2, 



(5) 



where i — (x + 8 x ,y + S y ,z), j — (x,y,z), and the displacement vector between spins i and j is 6 = (±1,0,0) or 
(0, ±1,0). In the following we will consider the forms of the driving field F which are perhaps the most relevant 
experimentally. Our main focus is the case of 'shear-like' linear variation of driving field with z, and the field acting in 
the x-direction only, such that the field components are F x (z) = 72, F y {z) = 0. Thus exchanges along x are enhanced 
or suppressed, while exchanges in the y and z directions proceed with equilibrium rates (AVF = 0). We also study 
the case of a V-shaped spatial dependence, F x (z) = 7 \z\, F y (z) = 0, such that the drive acts in the same direction 
throughout the system. Some results for a spatially uniform driving field in the x direction, F x = f = const, F y = 0, 
and a step-like field, F x (z) — sgn(z) • /, are also included. 

We have carried out extensive Monte Carlo (MC) simulations of the above model using single-spin and multi- 
spin 30, 31 coding techniques, in the latter case generalizing the driven multi-spin algorithm used previously to 3d 
systems. The multi-spin method we have adopted allows simulation of 64 independent systems (and hence greatly 



enhanced statistics) on a 64-bit computer system, using efficient bitwise operations. The state of an Ising spin may 
be represented by one bit; thus the values of a particular site in 64 different systems can be stored in a 64-bit variable 
(the natural word size). Bitwise operations operate on all bits of a variable at once, and are computationally cheap 
instructions. By combining these operations appropriately, and generating random bits with the required probabilities, 
the desired acceptance rates can be produced. 

Single-spin results provide a check as to the correctness of the multi-spin implementation. Parallelization via 
domain decomposition was employed to speed up simulations; the lattice was sub-divided along the x direction, and 
appropriate synchronization used when exchanging spins on and near the boundary between two domains. Data 
processing was also parallelized, owing to the large quantity of data produced by the multi-spin algorithm. Results 
reported here are from total run lengths of 2 — 4 x 10 7 MC sweeps (L x x L y x L z trial moves), the slow evolution of 
Kawasaki dynamics to a steady state proving to be less of a problem in 3d than in 2d [27] . 

The majority of the results shown here are for a system size L x = L y = 128 and L z = 10 or L z = 20, at a 
temperature T/T c = 0.75, where T c sw 4.5115 (f3 c = 0.2216544(3)) is the bulk critical temperature of the equilibrium 
3d Ising system [32] - We have checked that going to larger values of L x = L y < 192 has a minor effect on the results 
only for the smallest considered wall separation L z , i.e., as long as the longitudinal correlation length £n, which grows 
exponentially with L Zl is less than L x ,L y . We have also varied the temperature, firstly to investigate the effect of 
an increase to T/T c — 0.90, and secondly to study the behaviour near and below the roughening transition. For an 
equilibrium system in the thermodynamic limit, the roughening temperature is Tr « 2. 454 J |14j . In a finite system, 
the (pseudo-)transition will occur at a shifted value of temperature, which will be governed by the system size [33] . 
As temperature is increased in the smooth regime, either the interfacial width will reach the scale of L Z) or the lateral 
correlation length will reach L x or L y - in either case, the interface reaches the rough regime. We have thus covered 
a range of temperatures around the roughening temperature in the simulations, from T/T c — 0.4 to T/T c = 0.6. The 
roughening transition belongs to the universality class of the Kosterlitz-Thouless transition |34j . The renormalization 
group method of Kosterlitz 35 showed that the correlation length £n diverges very rapidly as T — ¥ Tr from below: 



£|l = Acxp 



B 



^/{Tr/T - 1) 



(6) 



where A and B are non-universal parameters. The numerical values for these constants obtained from MC simulation 
studies of the Ising interface in 3d are A = 0.80(1) and B = 1.01(1) [E]. The shift of the pseudo-roughening 
temperature can be estimated from the condition £|| ~ L x (= L y = L), which yields (Tr — Tr(L, L, L z ))/Tr(L, L, L z ) = 

AT ~ (B / \n(L / A)) . This is a very weak dependence on L and for our system size it gives AT ~ 0.04. At the same 
time the width of the interface diverges upon approaching the bulk roughening temperature, as w 2 ~ ln(£n). The 

condition w ~ L z yields AT ~ (B /(L 2 , — In A)) , which is a much stronger dependence on the finite dimension L z 
than that obtained from the condition involving £||. For L z = 10, AT ~ 10~ 4 , therefore we conclude that the shift of 
the roughening transition is governed by £m . Using that estimate of AT gives Tr(L, L, L z )/T c ~ 0.52, so the range of 
simulation temperatures should include the equilibrium pseudo-transition temperature. 

III. RESULTS 
A. Structure 

We first investigate the interfacial structure of the driven 3d Ising model in order to see whether also in three 
dimensions the effective action of drive is to increase the confinement of the interface. 

1. Magnetisation profiles 
The magnetisation profile along the z axis is calculated as 

LxL y \*,v I 

where the angles denote an average in the steady state. In a phase separated system m(z) changes sign across 
the interface, and attains values close to ±1 near the upper and lower walls respectively. In Fig. [2] we plot the 
magnetisation as a function of the scaled variable z = 2z/L z . Upon applying either shear-like or V-shaped drive, the 



magnetisation profile becomes 'sharper': m(z) changes sign more rapidly in the interfacial region, and there is a more 
extended flat region near either wall. The size of this effect increases with increasing 7. These trends are the same as 
in two dimensions [5S] , but for given driving strength, we find that the magnitude of the effect is smaller in 3d. 




FIG. 2: (Color online) Magnetisation profiles m(z) as a function of a scaled coordinate z between the walls. The system 
is L x = L y — 1 28, L z — 20, at a temperature T/T c = 0.75. Results for equilibrium simulations (7 = 0) with Kawasaki 
dynamics are shown (by symbols and lines), as well as for shear-like drive with several values of 7 (with symbols only). Results 
of rescaling to equilibrium are plotted in the inset; rescaling factors are ai = 0.88(5), 0.81(46), 0.7f (8), 0.61(4), 0.5f (28) for 
7 = 0.25, 0.5, 1, 2, 5, respectively. Error bars are of order or smaller than the line thickness or symbol size, and are not shown. 



It is possible to rescale the driven profiles to collapse back onto the equilibrium result: see Fig. [2] We interpret 
this as the drive acting to reduce the effective distance between the walls from L z to L*, and thus to effectively 
increase the confinement of the system. This increase may be quantified by introducing the scaling m(aj_z) « m eq (z), 
where aj_ = L*/L z is the ratio of the effective and actual wall separations. For shear- like drive with 7 = 1.0, we find 
a± = 0.71(8), whereas for two dimensions [3S], for the same value of 7 (7 is named ui there), a 2 ^ — 1/3.4 = 0.29. One 
may more formally express this in terms of a finite size scaling function for the profile. The scaling relation for an 
equilibrium fluid confined between two walls was proposed by Fisher and de Gennes [33] • In the absence of a bulk 
field it reads: 



m(z, T,L 



z/cq 



m h {T)M 



cq 



L, 



UTYUT) 



m b (T)M 



cq 



Lz'tbCn 



(8) 



M oq and M eq are 



where £(,(T) is the bulk correlation length, and m\>(T) is the spontaneous magnetisation in bulk, 
finite-size scaling functions: A4 Qq (u,w) is obtained from A4 cq (u,w) by changing the first scaling variable u = uw. 
Thus in equilibrium the shape of the scaling function can be varied by changing the wall separation at fixed T or, 
equivalently, by changing the temperature at fixed L z . We find that driving changes the shape of the interfacial profile 
at fixed temperature and L z such that 



m(z,T,L z ) 
m h (T) 



M 



cq 



L,* 



Lt'tbCn 



+ M co „(z) 



with L* < L, 



(9) 



where M. co „ is a boundary correction that decays away from the walls on the scale of £(,. Relation pi admits (in the 
scaling regime) another interpretation of the result fcty, namely, as the drive acting to reduce the effective temperature 
of the system at fixed actual distance between the walls. The estimation of the rescaling factor a± is obtained by 
rescaling the driven data to spline-interpolated equilibrium curves and minimizing the associated chi-squared statistic 
- the value of a± at the minimum is the optimal value. For the L z — 20 system, the 8 points in the centre of the system 
were included in the rescaling procedure - points further out are subject to stronger wall interaction effects. Both 
ways of driving yield very similar rescaling factors; more pronounced differences are observed for smaller values of 7 
and lower temperatures with the conclusion that effective confinement is stronger for the V-shaped drive. Rescaling 
fails for lower temperatures closer to and below the equilibrium bulk roughening transition temperature. 



The magnetisation profile may also been used to study the behaviour of the interfacial width. We measure the width 
via the second moment of dm/dz and study its variation with driving strength, wall separation L z , and temperature. 
Upon increasing driving strength 7 or /, the width reduces, as expected from the results for the full profile. For 
shear-like drive, we are also able to obtain data collapse for the behaviour of w/\/L z as a function of a scaling variable 
9 = L z ^/ S . Here s is an adjustable exponent. The division of the width by \JTT Z corresponds to the expected equilibrium 
behaviour [4j [15j [18], so that for 9 — > 0, w/ ' \f~L~ z — > constant. The scaling behaviour of the width is shown in Fig. pi 
for fixed temperature T/T c = 0.75, and a variety of wall separations and drive gradients in the ranges 10 < L z < 20, 
< 7 < 2. Previously we obtained data collapse as a function of the same scaling variable in the 2d system p6J (there 
the width was scaled by the 2d equilibrium behaviour, w ~ L z ). Remarkably, we obtain data collapse for the same 
value of exponent s = 0.3 as for the 2d system. From Fig. [3] we see that for small 7 at the larger values of L z = 16 
and 20, the data collapse is lost - we believe that this is because one moves out of the confined regime with w ~ \[L~ Z 
for these parameters. For these wall separations, the longitudinal correlation length £y becomes comparable to the 
linear dimension L of the interface, and the system crosses over to the regime where the dominant length scale is L. 



This does not require a large increase in L z , because £y ~ exp(«;L z /4), where the transverse length scale 



& 



for Gaussian interface fluctuations [TJ [TS]. Data collapse is regained for larger values of 7, because the effective wall 
separation L* < L z is the controlling length scale (from the discussion of the magnetisation profile above) , and L* z is 
small enough for the system to be in the "confined regime" . The inset of Fig. [3] shows the variation of the width with 
drive gradient 7 for shear-like and V-shaped drive. The trends are rather similar, with the width for given 7 very 
slightly smaller for V-shaped drive - this is consistent with the previous conclusion that the confinement is stronger 
for this drive type. 
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FIG. 3: (Color online) Scaling behaviour of the interfacial width w at a temperature T/T c = 0.75, as a function of the scaling 
variable 9 = L^ 8 , with s = 0.3. The width is scaled by y/L^ according to the equilibrium relation w/^/L^ = const., as 
discussed in the text. Different point types correspond to differing values of L z , from 10 to 20, as indicated. Inset: variation 
of the width with drive gradient 7, for shear-like drive (filled circles) and V-shaped drive (open circles). 



2. Spin-spin correlation functions 

In order to characterize the behaviour of the sytem on the two-body level, we first consider the microscopic interface 
pair correlation function, i.e., the spatial spin-spin correlation function at the midplane, defined as 

G(x, y, z = L z /2) = j±j- I J2 v{x', 2/', L z /2)a(x' + x,y' + y, L z /2) \ , (10) 

" \x ! ,y' I 

which depends on separations in both the x and y directions. Comparing the specific cases G(x, y — 0) and G(x = 0, y) 
provides information on the anisotropy in the x and y directions. In equilibrium, G(x 1 0) — G(0,y) for L x — L y , but 
for driven systems, the two functions may differ. Results in Fig. Plk show that shear-like drive causes both G(x,0) and 
G(0,y) to decay more quickly and for larger separations to saturate at larger asymptotic values than in equilibrium. 
In the x direction, this finding is in agreement with hydrodynamics results [23]; however in that study, the correlation 
length in the y direction was found to increase under shear, contrary to the trend in our system. We defer further 
exploration of this difference to the discussion of the height correlations below, since the height variable provides a 
more direct point of comparison between the systems. 

As in the 2d case, the spin correlation functions at intermediate separations may be transformed to the equilibrium 
result via a rescaling of the lateral coordinate, x or y: G(afiX, 0) « G cq (a;,0) and G(0, affy) ps G cq (0,y) [26]; see the 
inset of Fig. Eh for rescaling results for G(x,0). The rescaling factors are obtained via the same method as for the 
magnetisation profile; in this case, very small values of x or y are cut-off in the procedure, as are the tails of the 
functions, so that the rescaling procedure is carried out over 2 < x, y < 16. The a\\ parameters may be interpreted 
as ratios of lateral interfacial correlation lengths in and out of equilibrium: aif = £jf /^' cq . Rescaling the driven 
results produces a?, < 1 and ajf < 1, so that the correlation length is reduced under drive - this tallies with the 
faster decay evident from Fig. Ha, which shows results for equilibrium and shear-like drive. The x-y anisotropy may 
be measured by the ratio aYJaf,; this is consistently slightly smaller than unity, leading to the surprising conclusion 
that the correlations are slightly more suppressed in the y (vorticity) direction. As with the magnetisation profile, 
the effect of drive is much weaker in three dimensions than in two - for example, in 2d, ajf = 1/1.27 for 7 = 0.025, 
while in 3d, aff = 1/1.271 ss 0.79 for 7 — 0.25: a ten times larger field gradient is required to produce a comparable 
confinement. As a result, the rescaling procedure works for much larger values of 7 than in 2d - a stronger drive is 
required to push the system into a regime which is too far from equilibrium for rescaling to be possible. The V-shaped 
drive has a similar effect on the interfacial correlations: for example, with 7 = 0.25, a?, = 0.74, corresponding to a 
slightly stronger confinement effect than with shear, consistent with the findings for the magnetisation profile. 

Fig. [4]j shows the effect of varying the temperature on the interfacial correlations. Lowering the temperature to 
T/T c = 0.5 (below the roughening transition in a bulk equilibrium system) in equilibrium results in correlations 
G(x, 0) that decay only to ~ 0.6 for the largest separations. Since on the lattice the average interface position 
lies between two lattice points (for zero overall magnetisation), one measures the correlations just either side of the 
interface (which side does not matter, due to symmetry). Thus at zero temperature, G(x, y,z — L z /2) — 1 for all x, y, 
since the interface is perfectly flat at T = 0. This explains the observed increase of asymptotic values of G(x, 0) for 
low T. Moreover, the width of smooth interfaces is of order of the bulk correlation length, which at low temperatures 
is ~ 2-3 lattice spacings. Therefore, for low temperatures G(x, 0) essentially measures correlations in the bulk-like 
phase. We also see from Fig. [4b that driving the system enhances the asymptotic value further for T/T c — 0.5, which 
indicates that at fixed temperature the bulk-like phase is more ordered under drive than in equilibrium. As for the 
magnetisation profile, rescaling does not work for the low temperatures - the drive affects the asymptotic value more 
than the decay rate for these temperatures. 

3. Height-height correlation functions 

We now turn to the interfacial height-height pair correlation function. The interface height is obtained via a 
coarse-graining procedure, which produces a single-valued height function h(x,y) from the real microscopic configu- 
ration. The latter contains 'bubbles' of one phase in the other, and 'overhangs' at the interface, meaning a height 
cannot be defined from it directly. To coarse-grain, we use a simple method which we found to be successful in 
the 2d system, where it gives results equivalent [26j to a more complicated coarse-graining method |37j . The height 
h(x,y,t) = (l/2)^2 z cr(x,y,z,t) is simply a sum of the spins over a column. Thus when there are equal numbers 
of '+' and '-' spins in a column, h(x, y) = 0, while when there is a majority of one species, h 7^ 0. The general 
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FIG. 4: (Color online) Spin-spin correlation functions G(x, y = 0) and G(x = 0, y) as a function of scaled separation x/L x or 
y/Ly. (a) Results for a 128 x 128 x 10 system at T/T c = 0.75. The equilibrium result is shown, as well as driven results for 
shear-like drive at several values of drive gradient 7, as indicated. In the inset the driven results for G(x, y = 0) are rescaled via 
the parameter a|, as described in the text; rescaling factors are af = 0.95, 0.79, 0.69 for 7 = 0.05,0.25,0.5, respectively. The 
maximum separation for the correlation functions is x/L x — y/L y = 0.5, due to the periodic boundary conditions; beyond the 
displayed region the functions have essentially reached their asymptotic values, (b) Results for G(x, y = 0) for 128 x 128 x 10 
systems at T/T c = 0.75 and T/T c = 0.5, for equilibrium and two strengths of shear-like drive. Below the equilibrium roughening 
temperature, i.e., for T/T c — 0.5, the asymptotic value is large, and increases with stronger drive. 



height-height correlation function depends on spatial separations x, y and on temporal displacement t: 

C(x, y, t) = -^- ( Y, K*', V', t')h(x' + X, y'+y,t' + t)), (11) 



\x',y' 



where the angles indicate an average over time. We first consider the equal-time correlations, with one of the spatial 
separations set to zero: Fig. [5] shows results for C(x,y — 0,t = 0). In two dimensions, C(x,t = 0) (also G(x)) in 
equilibrium exhibits strong anti-correlated regions for medium-to-large separations, presumed to be finite size effects 
[27j . These are not present in 3d for the system sizes considered - the functions decay to zero without becoming 
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signifcantly negative, indicating less severe finite size effects; an explanation may be the following. With conservative 
dynamics, a positive-height 'bump' must be accompanied by one with negative height, since Y] h = 0. In 2d, these 
must lie on the same x-z layer (the only one), so an anti-correlation is measured in C(x). However in 3d there are 
L y x-z layers, so the pair may be located in different layers, meaning C(x,y — 0) does not necessarily display anti- 
correlations. Turning to the driven cases, we see that applying shear-like drive leads to more rapid decay of C{x, 0, 0), 
as well as a smaller initial value C(x — 0) (a measure of the interfacial width). The magnitude of this effect increases 
with increasing 7. In this section, the results shown are for shear-like drive, but the conclusions also apply to V-shaped 
drive - as for other static quantities, the effect is similar to shear, with slightly greater correlation-suppression. 

As for other quantities, rescaling to the equilibrium result is possible in the same manner as in 2d. For the height 
correlations, the rescaling takes the form a]_ C{af<x) ss C cq (x). The values of a± and a\\ are those obtained from 
the rescaling of the magnetisation profile and the spin-spin correlation function for given simulation parameters. In 
2d this procedure was motivated by Weeks scaling in equilibrium [3]: C eq (x) « w 2 C{x/^, q ), where C is the scaling 

function, and since w oc L z in 2d, the correct rescaling involves a\. Although Weeks scaling does not hold in 3d, this 
procedure works reasonably well for 7 < 0.5, except at zero separation, where the rescaled values become larger than 
the equilibrium width. As with the spin correlations, this range of validity is much greater than in two dimensions. 
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FIG. 5: (Color online) Height-height correlation function C(x, y — 0) as a function of separation x, for a 128 x 128 x 10 system 
at T/T c = 0.75. In the inset, data for non-zero drive are rescaled to the equilibrium result via the relation aJ 2 C{a»x) w C eq (a:) 
given in the text, where the values of a± and ojf are obtained from the rescaling of the magnetisation profile and spin-spin 
correlation functions, respectively. 



Furthermore, we are able to fit the results for C(a;,0,0) and C(0,y, 0) for small-to-intermediate separations to the 
equilibrium capillary wave result for the height correlation function in 3d; see Fig. [6] for results for shear-like drive. 
This procedure was used previously by Derks et al. to describe their experimental data, where an excellent fit was 
obtained [2U. By integrating Eqn. S, the (equilibrium) capillary wave result for C(x, 0, 0) in d — 3 and in the limit 



L x — _L„ ^ 



is ® 



C(x,0,0) 
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(12) 



where Kq is the modified Bcssel function of the second kind. The upper wave-number cutoff has been sent to infinity 
in order to obtain an analytic result; in order to regularize the integral at x = 0, a shift A is introduced. Bedeaux 
and Weeks j3| give A ss l/(<7max£ii), where g max is the original wave-number cutoff in the integral ([2]). Combining (12) 



with the capillary-wave result for the interfacial width, 



w 2 = C(0,0,0) 
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(13) 



we are able to substitute for the (unknown) q max in terms of the width, and obtain a fitting form with two parameters: 
the correlation length £h and the pre-factor h^T/V. In the above, we have specialised to separations in x rather than 
the radial distance r usual in CWT, since isotropy is broken in the non-equilibrium situation. For separations in y, 
the form for C(0,y, 0) is the same, but different values of the parameters are expected - i.e., the correlation length 
(£?l ) will be different, as will the pre-factor. The intepretation of the latter quantity is difficult. Indeed, the interface 
tension is an equilibrium concept and cannot be carried over directly to non-equilibrium situations, so the meaning 
of the pre-factor is not initially clear - here we just note its anisotropy. 

The equilibrium fit wanders off the data for larger separations; this may be due to a (less serious) manifestation of 
the finite size effects encountered in 2d, which were mentioned above, and the conserved order parameter. (This is 
most obvious for the equilibrium data on the log scale in Fig. [6j where the data diverge as they approach zero and 
become negative). For the driven cases, as drive becomes stronger, the fit works for a smaller range of separations - 
the example of shear-like drive is given in Fig. [6] This trend is expected from the findings for the rescaling procedures 
applied above - initially the system is "close enough" to equilibrium for CWT to be approximately applicable, but 
as drive increases, this ceases to be true. From the fits we obtain the equilibrium and non-equilibrium correlation 
lengths in the x and y directions, £jf and £?,- see the inset of Fig. VM for their variation with 7. The trend of decreasing 
correlation length with increasing drive strength mirrors the one found in 2d, although there we were not able to 
obtain the correlation length reliably, due to the difficulty of fitting the correlation function data over a reasonable 
range. We also note that £}! is consistently smaller than ^rF, in agreement with the earlier conclusions, based on the 
behaviour of the spin-spin correlation functions, that correlations are slightly more strongly suppressed in the vorticity 
direction than in the driving direction. 
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FIG. 6: (Color online) Fits of the equilibrium and non-equilibrium height correlation data for C(x, y = 0, t — 0) to the 
asymptotic capillary wave prediction given in the text. The system parameters are the same as in Fig. [5] Inset: correlation 
lengths £« and ff, (via fitting C(x = 0, y, t — 0) data) obtained from the fits, as a function of shear-like drive gradient 7. 



The suppression of correlations we find in the drive (x) direction is in agreement with the hydrodynamics work of 
Thiebaud and Bickel [25], who studied phase separated fluids between two walls under shear. This trend is also the 
same as in the 2d Ising system, and both microscopic (G(x), spin-spin) and coarse-grained (C(x), height) measures 
of correlations give the same conclusion. Both the theoretical and simulation findings disagree, however, with the 
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experimental results of Derks et al. 24 , who found an increase of correlation length in the flow direction when shear 
was applied to a phase-separated colloid-polymer mixture. The fact that we have used the same method of fitting 
the height correlation data to the equilibrium capillary wave form as Ref. [53] , makes the method of comparison the 
same, at least. 



Finally, we consider the pre-factor resulting from the fit of the height correlation data to the CWT form (12 1. In 
equilibrium, the pre-factor is proportional to ksT/T, where T is the surface stiffness, which for a continuum fluid is 
the interfacial tension the free energy associated with the interface. Out of equilibrium, this quantity is not defined, 
so the meaning of the pre-factor resulting from the fit is not clear. Numerically, we find that the pre-factor from 
fitting C(x, 0, 0) is a decreasing function of drive gradient 7 for shear-like drive. If one defines a "non-equilibrium 
surface tension" from the CWT fit, and further takes the temperature T to be fixed (i.e., the value of the parameter in 
simulations), the conclusion is then that this "tension" increases as the system is more strongly driven. This procedure 
of defining an effective non-equilibrium surface tension via the CWT fit was the approach taken in the analysis of the 
experiments of Derks et al. 24J, who also found this tension to be an increasing function of shear rate. The CWT 
fits of experimental data in Ref. 24J show that an increase in a "non-equilibrium surface tension" is accompanied by 
an increase in correlation length, but our simulations show the opposite relationship - a decreasing correlation length 
as the effective interfacial tension increases. Our findings seem to be inconsistent with the CWT result (see Eqn. (pi) 
and below) £11 oc vT. However, in our system at equilibrium £11 oc vTexp(L z /(4£f,)), so that the effective increase of 
the confinement due to driving (reduction of L z ) wins over the effective increase of T. 

4- Structure factor 

To further complicate the situation, our finding of a decrease of correlation length in the vorticity (y) direction is 
in disagreement with Ref. |25j . where an increase was found. Intruigingly, Thicbaud and Bickel found the structure 
factor S*(q) = S(q x ,q y ) to be unaffected in the q y direction by the application of shear |25j . The same was concluded 
for the uniformly driven system from the analytic approach based on the time-dependent Landau-Ginzburg functional 
by Leung [22] . The static structure factor in our system is accessible via a two-dimensional spatial Fourier transform 
of the equal-time height correlations, C(x, y, t = 0): 

S(<i) = \F{C(x,y 7 t = 0)}\\ (14) 

where T {} denotes a two-dimensional spatial Fourier transform, q = (2n x Tr/L x ,2n y ir/L y ), and n x , y = 
0, 1 . . . ((L XtV /2) — 1), so that q x and q y lie on the range . . . (ir — (2n/ 'L x ^ y )). In Fig. p] we plot l/S{q) for equi- 
librium and driven systems, along either the q x or q y direction, as a function of q 2 = |q| 2 . From Eqn. pi), in 

equilibrium, one expects l/S(q) oc (T/k B T) (q 2 + £f 2 )- In Fig. 7 we fit the equilibrium data to this form. The data 
shown are along the direction with q y = 0, although we have checked that the equilibrium structure factor behaves 
the same along the q y direction, as expected. For q x < 1, this behaviour is indeed observable in the simulations, but 
f° r Qx ^ 1, the data diverges from the CWT prediction, when other powers of q x presumably become important. 

Turning to the non-equilbrium behaviour, we see that for shear-like drive, S(q) is affected (suppressed) in both the 
drive (x, blue crosses) and vorticity (y, red squares) directions, but the effect is smaller in the vorticity direction. 
These results are in disagreement with the hydrodynamics results [25] ; however, since the effect in the drive direction 
is stronger than in the vorticity direction, the latter effect could possibly be of higher order than was considered in 
Ref. [25]. The data for shear- like drive along q y = are also fit to the equilibrium CWT form in Fig. pi we see 
that as for equilibrium, the fit is reasonable for q x < 1. Additionally, the intercept at q x = is greater, indicating 
a smaller lateral correlation length, as found in real space above. Indeed, one can compare the parameters resulting 
from the CWT fits in real space and Fourier space. We find that the qualitative trend for the correlation length is the 
same, but do not obtain quantitative agreement - the values obtained from the real space fits are consistently larger. 
These differences are expected - for equilibrium, they can be caused by the finite system size and lattice discretization 
effects. For non-zero drive, the effect of deviations from CWT can be different in real and Fourier space. Additionally, 
the fits in Fourier space are for small q (long wavelengths), while the real-space fits are for small separations, so the 
length scales the fits are applicable to is not necessarily the same. For V-shaped drive, we find that for given drive 
gradient 7, the results are similar to those for shear-like drive, with slightly greater suppression of the structure factor 
at small q. 
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FIG. 7: (Color online) Structure factor S(q), as defined in (141. Data are shown for S(q x ,q y = 0) and S(q x = 0,q y ), as a 



function of the squared norm of the wave vector, q = q x or q y , respectively. As before, the system dimensions are 128 x 128 x 10, 
and the temperature is T/T c — 0.75. Equilibrium data are shown, as well as for shear-like drive with 7=1, and fits to the 
CWT form are displayed for equilibrium and for S(q x , 0) with 7 = 1, up to q x = 1. 



B. Capillary wave transport 



We now consider the dynamics of the height-correlation function defined in (11). The primary limit of interest is 
C(x, y = 0, t) which we find shows evidence of capillary wave transport along the driving direction, for suitable forms 
of the current profile. Previously (in the context of d = 2) [3S], we conjectured that the capillary wave fluctuations 
on an interface will be transported by an external driving field, provided that the lateral order parameter current 
has a component which is an odd function of distance from the interface. Thus only purely even current profiles are 
expected to give no transport. 



1. Order parameter current 

In the 3d system, the order parameter current profile is defined as j(z) = j+(z) — j__ (z), where j a (z) is the current 
profile of the ±1 spin species. The components j%{z) and j%(z) oi j a (z) are the net number of spins of that species 
moving in positive x and y directions per unit time at perpendicular coordinate z. We also note that the Ising 
symmetry means that the order parameter current can also be written as j(z) = 2j + (z); the first definition may be 
applied to systems lacking the Ising symmetry, i.e., liquid-gas or liquid-liquid interfaces. 

As shown in Fig.^1 the shear-like drive F x (z) = jz, F y (z) — gives a purely odd order parameter current component 
j x (z). Since the y component of the field is zero, j v (z) = 0. For the V-shaped drive, the current is an even function 
of z, since F x (z) = 7 \z\. We thus expect transport along x for the shear-like drive, but none for the V-shaped drive. 
For the same value of 7, the currents for the two drive types almost coincide in the region z > 0, where the driving 
fields are the same in magnitude and direction. 

In Figs. 8(a) and |8(b)| j x (z) is shown for various drive gradients 7, for temperatures above and below the (equilib- 
rium) roughening transition. Looking at the high temperature data in Fig. 8(a) we see that for small 7, |j x (5;)| has 



maxima at the walls; as 7 is increased, plateaus develop with the maximum current shifted slightly from the wall. 
Eventually for strong drive the maxima become localized near the interface. This reflects the competing effects of local 
drive strength and current carrier availability (H — pairs): for large 7, the drive strength is essentially saturated at the 
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walls, so the greater carrier density at the interface eventually becomes more important. Below the bulk roughening 
temperature (T = 2.4, Fig. 8(b)| , the current |j :r (^)| also has maxima at the walls for small 7, and quickly develops 
maxima at the middle two layers as 7 is increased. These maxima appear for much weaker drive (approximately six 
times smaller 7) than they do for T/T c = 0.75. They are also localised to the two middle layers either side of the 
interface, and are much more pronounced than at the higher temperature; this indicates that at low temperatures 
the interface region is very sharp, reduced to approximately two lattice spacings. For strong drive, the greater carrier 
density at the interface again 'wins', and these maxima become global. We also note that Ij^z)! is roughly five times 
smaller than that at the higher temperature, since the carrier density is much smaller due to the increased bulk and 
interfacial order. 



Finally, Fig. 8(a) also shows an example of mixed symmetry in the current profile. The driving field is of the 'V 
type, but with different values of 7 in the upper and lower halves of the system: 7/ — 0.25 in the lower half, j u = 0.5 in 
the upper. Thus the total driving field can be written in the form F{z) = 71 \z\ + j 2 z, with 71 = 0.75/2, 72 = 0.25/2, 
showing the even and odd components explicitly. The current profile reflects the asymmetry in the drive: in the lower 
half of the system, the current matches that for a (symmetr ic) V -shaped drive with 7 = 0.25, while in the upper half, 



it matches that for either V or shear with 7 = 0.5 (see Fig. 8(a) ) - the crossover occurs over a single lattice spacing 



2. Space-time correlations. 

We investigate whether the conjecture for the occurrence of capillary wave motion holds for 3d systems by measuring 
C(x, y = 0, t) for different forms of driving field, which produce differing current profiles. For current profiles with 
odd symmetry in z (or more generally, profiles with an odd component) , we expect to see evidence of capillary wave 
transport in C(x, y = 0, t). Fig. M shows that this is indeed the case - see the main panel for results for C(x, y = 0, t) 
for shear-like drive, and the inset for V-shaped drive. For time difference t = 0, the peak lies symmetrically around 
x = due to the translational invariance ensured by the periodic boundaries along x. However, at time differences 
t > 0, the peak moves to negative x values for shear-like drive, indicating that now the greatest correlations are 
between spatially-displaced points. We interpret this to mean that wave-like height fluctuations are being coherently 
transported along the interface by the drive. For the V-shaped drive, the peak remains at x = for all t, showing the 
absence of wave motion. In both cases, correlations decay with increasing time difference, due to thermal noise. We 
note that the rate of decay of correlations is much faster for driven systems than it is for equilibrium systems with 
Kawasaki dynamics. 

We have also investigated other forms of driving field, for example spatially uniform driving field in the x direction, 
F x = f — const, F y = and step-like drive: F x {z) = sgn(z) • /. The former produces an even order parameter current 
profile whereas the latter an odd one, and results for C(x, 0, t) (not displayed) show that wave motion does not occur 
for a uniform drive but occurs for a step-like drive, consistent with the conjecture. For the case of the asymmetric 
V-like drive discussed in the previous section, which has mixed symmetry, we expect to see wave movement, since the 
current profile is not purely even, but like the driving field itself, can be written as a sum of even and odd components. 
Indeed we find this is the case, with the peak of C(x, 0, t) moving with time. From these results we conclude that the 
criterea for capillary wave motion are the same in the 2d and 3d Ising systems. 

Having established the occurence of wave motion, it is natural to investigate the dependence of the wave velocity 
on system parameters. We measure the speed of the peak of C(x,0,t), u pC ak, and vary the driving strength (7 for 
shear- like drive, / for step-like drive), temperature and wall separation L z . As shown in Fig. |10[ w n nak shows linear 



^pcak 



variation with 7 or / for fixed temperature and system size, for 7, / < 2. For shear-like drive, the gradient of v pea k(7) 
is close to 2 in this range. We also see that varying L z has a rather small effect on the peak velocity - doubling 
L z from 10 to 20 reduces the gradient of w pea k(7) by only a few percent. Changing the temperature from 0.75T C to 
0.90T C also has small effect, in the other direction - the peak moves faster for the higher temperature at given L z , 7. 
For the step drive, w pe ak also seems to be linear in driving strength / for small /, with a reduced gradient compared 
to shear-like drive. For both forms of drive, non-linearity appears to set in for 7 > 2. For the mixed symmetry case, 
where the driving field could be written as F{z) = 71 \z\ + 72Z, we find that the velocity of motion is smaller than 
that for a purely odd field with 7 = 72 - the velocity in the mixed case is linear in 72 and approximately 80% that 
of the pure case for 7, 72 < 2. For lower temperatures near and below the equilibrium roughening temperature, we 
find that the interfacial motion still occurs, with a much reduced velocity; correlations also decay much more quickly 
with time. 
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3. Dispersion relation 

In order to characterize the dynamics of capillary waves we consider the evolution of the spatial Fourier modes 
h(k,t) of the height function h(r,t) :— h(x,y,i) defined by: 

L 

Mr,i) = ]Te 2 ™( k / L > r Mk,i), (15) 

k=l 

where k = (k x , k y ) and r = (x, y) are two-dimensional vectors and the division k/L is defined as k/L = (k x /L x , k y /L y ). 
(Recall that x and y are integer coordinates of the lattice.) The sum denotes summations over integer k x and k y from 
1 to L x and L y , respectively. Because /i(x,y,£) is a real function, the Fourier transform has the conjugate symmetry 



h(k x ,k y ,t) = h*{L x — k x ,L y — k y ,t), i.e., there are L x /2 x L y /2 independent terms in (15 1. Each complex Fourier 
component /i(k, £) can be written in terms of its modulus and phase 0(k, t): 

fe(k,t) = |fc(M)|e* (k,t) , (16) 

where — tt < (f>(k,t) < tt. If 0(k, t) = uj(k x ,k y )t — [2iik x / L x )v x (k y )t, each mode would correspond to a travelling 
wave moving in the x direction with a velocity v x and h(x, y, t) = h(x + v x t, y) (with no dispersion). In a steady state, 
the phase shift of each mode is a fluctuating quantity, with a distribution which when measured with increasing time 
intervals spreads and decays quickly to zero. However, at short times, we are able to measure its mean value in unit time 
to obtain the dispersion relation of the frequency cj as a function of the wave vector q = (q x , q y ) = (2itk x /L xl 2irk y /L y ) 
as: 

«(q) = aig(h*(ci,t)h(q,t + dt))/dt. (17) 



Results for w(q) = uj(q Xl q y ) calculated for dt equal to 1/10 of an MC sweep are shown in Fig. 11 as a function of q x 
for several values of q y . We observe similar behaviour for all considered temperatures, including below the roughening 
temperature (data not shown). The shape of uj(q x ,0) is very similar to that obtained in the two-dimensional system 
(see Ref. [Ml), which suggests to use the same analytical formula for describing the dispersion relation: 

^iqxiqy — const) = (v + 2u)sin(q x ) — usm(2q x ) + ssin 2 (q x ) (18) 

with small-g^ expansion 

u(qx,q y = const) = vq x + sq x + (u- v/6)q x . (19) 



Indeed, as can be seen in the Fig. 11 our data for w(q x ,0) fit well with v = 0.0089(8), u = 0.01095(15), s = 0.0071(9). 



Waves with larger q y are less dispersive in the sense that the corresponding coefficients u and s are smaller. The 



fastest mode is the one with q y — n, or a wavelength of two lattice spacings; uj(q x ,Tr) can also be fitted using (18) 
with v = 0.0547(5), u = 0.00530(26), s = 0.0014(15). As we pointed out in Ref. [2SJ, if the dynamics of the height 
function are modelled by the following linear transport operator (propagator) 

L = (d t -vd x + udl), (20) 

then the Ansatz in the form of the travelling wave Aexp(i(ut + q • r)) yields the first and the second terms in 



the fit function (18). In (20), time is treated as a continuous variable whereas spatial derivatives are discrete: 
d x h(x, y, t) — [h{x + 1, y, t) — h(x — 1, y, t)] /2, and a 5-point stencil is used for d x . The third term, which for small q x 
gives the quadratic dependence in ( 19 ), can be obtained if one allows an imaginary contribution is(d x + d x /A) (with 3- 



and 5- point stencils) to the linear operator L. One may be able to understand the presence of complex coefficients in 
the transport equation for the height function h(r,t), by recognising that the plane wave solution above neglects the 
dependence of the amplitude A on the wave number. In fact the average modulus |/i(q, t)\ of each complex Fourier 
component varies significantly with q = |q|, even in the absence of driving (see the plot for the structure factor Fig.[7|. 
Taking this into account, it should be possible to derive an equation for the amplitude as well as the phase - this 
may aid in understanding the presence of real and imaginary parts in the transport operator. The t — > oo limit of 
the solution of the amplitude equation should yield the static structure factor, which may be compared to simulation 
data. We leave the interesting and difficult problem of deriving the full transport equation from these considerations 
to future work. 
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IV. CONCLUSIONS AND OUTLOOK 

We have presented evidence from Monte Carlo simulations that the main characteristics of intcrfacial structure and 
dynamics in the laterally driven stochastic lattice gas model are generic and persist to three dimensions. Far above 
the equilibrium roughening temperature the structure of the interface confined between two walls is affected by lateral 
driving in a way similar to that resulting from an increase of confinement (by reducing the distance between the walls) 
of the equilibrium system. However, the effect of drive is much weaker in three dimensions than in two, plausible 
due to the different nature of low-energy interfacial fluctuations, i.e., the "spike" -like excitations rather then interface 
"wandering" . In the case of the shear-like drive, our findings for the decay of the height-height correlation functions 
and for the structure factor are in partial agreement with recent results from fluctuating hydrodynamics |25j . The 
discrepancy concerns the behaviour in the vorticity direction. We have found a decrease of correlation length in this 
direction whereas hydrodynamic calculations predict an increase. Also, our data show that the structure factor is 
suppressed both in the drive and the vorticity directions. In Ref. |25| it is concluded that S(q x , q y ) is unaffected in the 
vorticity direction. Moreover, our results for the interfacial width are in agreement with the experiment of Ref. (24) . 
but the results for the interfacial correlation length in the flow direction are not. It would certainly be desirable to 
carry out more studies, experimental, theoretical and simulation-based, to clarify these discrepancies. Nearer (and 
also below) the bulk roughening temperature, the confinement/ordering effect of the drive is reduced, since there are 
no large fluctuations to "smoothen out" . The picture of an equilibrium system under a greater effective confinement 
no longer seems to apply for these low temperatures. 

The conjecture made in Ref. [35J for the occurence of lateral transport of interfacial fluctuations is supported by 
our results in three dimensions. The transport also occurs for low temperatures below the equlibrium roughening 
transition. However, the dynamics of the lateral propagation of capillary waves which gives rise to the observed 
dispersion relation is still not understood. A full treatment should include non-linear effects, and treat the time and 
wave-vector dependence of the Fourier amplitude. Further insight, especially into the effect of coupling between bulk 
and interfacial degrees of freedom, could be gained by considering theoretical models based on the order parameter 

EHJ. 
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FIG. 8: (Color online) Order parameter current profile component j x (z), for the system parameters L x = 128, L y = 128, 
L z — 20. (a) Temperature T/T c = 0.75. Results are shown for equilibrium (zero current), shear-like and V-shaped drive, and 
the case of mixed symmetry in the driving field. In the latter case the lower and upper-half 7 values are 7; = 0.25, j u = 0.5. 
For z > 0, the currents resulting from shear-like and V-shaped drive with 7 = 3 coincide, as do those from mixed symmetry 
and shear with 7 = 0.5. (b) T — 2.4. Results are shown for equilibrium, as well as shear-like and V-shaped drive with various 
values of 7. 
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FIG. 9: (Color online) Time-displaced height-height correlation function C(x, y — 0,t) for a 128 x 128 x 10 system at T/T c = 0.75. 
Main panel: snapshots at time displacements t, as indicated, for shear-like drive with gradient 7 = 0.5. The peak of the 
correlation function moves to the left with increasing time difference, indicating capillary wave motion. Inset: results for 
V-shaped drive with strength 7 = 0.5. Correlations decay without movement of the peak. 
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FIG. 10: (Color online) Speed of movement of the peak of C(x, 0, £), w pe ak, as a function of driving strength, 7 or /, for shear- 
and step-like drives, respectively. Shown is the effect of varying L z for fixed driving strength and temperature, and also the 
effect of increasing the temperature from T/T c = 0.75 to 0.9. 
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FIG. 11: (Color online) Dispersion relation u}(q x , q y ) for shear-like drive as a function of the wave number q x for several values 



of q y . Data are for the system parameters L x 
analytical formula (18 1. 



128, L y = 128, L z = 10 at temperature T/T c = 0.9. Solid lines are fits to the 



